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I .  INTRODUCTION 


The  Discrete  Fourier  Transform  has  traditionally  been  the 
dominant  tool  in  digital  signal  processing  for  extracting 
waveform  features  required  in  the  analyses  of  signals. 
Characteristic  of  the  Fourier  transform  is  the  global  manner 
in  which  it  operates  on  the  input  signal,  thereby  rendering 
this  technique  incapable  of  providing  time  localization  of 
frequency  components.  The  analyses  of  signals  that  are 
transitory  in  nature  will  particularly  suffer  the  deleterious 
effects  of  this  limitation  and  the  utility  of  the  discrete 
Fourier  transform  for  transient  signal  analyses  may  therefore 
be  diminished.  Various  data  windows  have  been  used  in  the 
past  in  an  attempt  to  improve  performance,  but  all  involve 
tradeoffs  in  resolution  and  sidelobe  levels  that  may 
unacceptably  affect  the  analysis.  [Ref.  l:pp.  63-82] 

Recent  appearance  in  the  literature  of  several  articles 
describing  wavelet  transforms  and  the  related  multiresolution 
analyses  have  created  interest  in  determining  their 
suitability  for  signal  analysis  applications;  specifically,  to 
overcome  the  time  localization  limitation  of  the  Fourier 
transform.  This  thesis  examines  the  potential  of  the  wavelet 
transform  for  signal  analysis  purposes. 
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II.  REVIEW  OF  WAVELET  THEORY 

A.  DISCRETE  WAVELET  TRANSFORM 

Wavelets  are  families  of  functions, 

^'*4*  (2-1) 

generated  by  the  dilations  and  translations  of  a  single 
function.  For  digital  signal  processing  applications  it  is 
necessary  to  discretize  the  parameter  values  and  fix  the 
initial  dilation  step  and  the  translation  step  bg  such  that 

ao>l,  bQ*o 

Equation  2.1  becomes, 

ba,„ix)  =af^h{aoX-iibo)  m.neZ  (2.2) 

with  3=30®  and  b^nb^ao'” 

The  translation  parameter  is  therefore  dependent  on  the 
dilation  parameter.  A  large,  positive  m  will  contract  the 
function  p  and  cause  the  translation  step  to  shorten  while 
a  large  negative  m  dilates  the  function  p  and  lengthens  the 
translation  step  size.  The  inverse  relationship  ensures 
coverage  of  the  entire  range.  [Ref.  2:  pp.  909-910] 
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Associated  with  the  discrete  wavelet  is  the  discrete 


wavelet  transform,  T,  which  maps  the  function  /  to  a  sequence 

indexed  by  Z^,  ^  ^  - - 

( Tf)  * <a?x-njbo)  f{x)dbc  (2.3) 


If 


[A(5)]M  <  (2.4) 


and  h  has  sufficient  decay,  and  if  T  has  a  bounded  inverse  on 
its  range,  the  set  <h^>  forms  a  frame  for  all  square- 
integrable,  one-dimensional  functions  f, 

f(x)€L^(m  (2.5) 

The  significance  of  a  frame  is  that  numerically  stable 
algorithms  exist  which  allow  f  to  be  reconstructed  from  the 
wavelet  coefficients  <h^,f>.  [Ref.  2:  p.  911] 

Selection  of  appropriate  h^a^,  and  is  influenced  by 
various  factors.  The  desire  to  take  advantage  of  previously 
published  work  and  to  minimize  redundancy  in  the  wavelet 
representation  resulted  in  choices  of  =  2,  bp  =  1,  and 
several  h  (discussect  below)  such  that  the  constitute  an 
orthonormal  basis  of  compact  support. 


B.  NULTIRESOLUTION  ANALYSIS 


Multiresolution  analysis,  as  the  name  implies,  describes 
the  function  f  as  a  series  of  approximations  of  the 
function  at  different  levels  of  resolution  and  consists  of  a 
family  of  embedded  closed  subspaces 

V„<=  (St)  t  ...  cV.2cVliCV’oCV'iCV'2 ...  (2.6) 

such  that 

and 

f(-)  (2.8) 

Also,  there  is  a  scaling  function  <|>(x)  ev^  such  that  the 

(where  <|>^^(x)  =2®''^<t)  (2®x-n) )  constitute  a  basis  for 

=  spaHT^^  (2.9) 

Let  P„f  be  the  orthogonal  projection  of  f  onto  V^.  From  the 
preceding,  it  can  be  seen  that  liin„.,Pj„f  =  f,  for  all  f€L^(R)  . 

[Ref.  3:pp.  967-968] 
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Considering  our  parameter  choices  and  the  decision  to  use 


orthonormal  bases,  if  we  define 

C^^f)  =  <<l)an,f> 


(2.10) 


then 

^  Vn,  (2.11) 

Now  define  =  P^if  -  and  as  the  orthogonal  complement 
of  (W„xV„)  so  that  =  V„©f4.  Then  is  the 

orthogonal  projection  of  any  f  e  (R)  onto  the 

orthogonal  complement  of  in  The  are  scaled 

versions  of  P/p  , 

f(-)  e  -  f{2«  •)  6  (2.12) 

and  the  are  orthogonal  spaces  which  sum  to  L-(R) 

lMR)=®P4  (2.13) 

Similar  to  the  there  exists  in  a  wavelet  function  ijr 

15*  E  • 

such  that 


5 


W„  =  span  I 


(2.14) 


where  ilraaj(Jc)  =  (2®x-n)  .  If  we  define 

d^(f)  (2.15) 


then 

=  <2.16) 

Naturally,  the  function  may  be  reconstructed  from  its 
projections  onto  the  bases  vector  spaces.  [Ref.  2;pp.  916-918] 
1.  Properties  of  the  Wavelet-Scaling  Function  Pair 

Following  is  a  brief  summary  of  the  relationship 

between  the  <t>  cind  i|r  families.  Since  (|>(,  (x)  eV'^cVi  ,  it  can 

be  written  as 

♦.(x)  =  j  (♦,(x),<ti.U-2))(f,<x-|)  (2-17) 

n—ei  ^  ^ 
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Let 


h{n)  =  J‘^<l>o(x)<|>o(2x-n)  dx 
=  2'^(4>,(x),4),(x-|)) 


so  that 


4»oU)  »  v/25;  ij(n)4)i(x-^) 


=  2  A(/I)<{>o(2X-/2) 


The  Fourier  transform  of  this  equation  is 

®  (2&))  =  Midi)  ®  (w) 


where  we  have  defined 

•• 

H{d>)  =  52  A(n)e-J"" 

n*-«* 


To  satisfy  the  above  relationships,  the  following 
must  hold: 

WiO)\  =  1 

|;f(u)|2  +  lff(a)+m)l2  =  1 


(2.18) 


(2.19) 


(2.20) 


(2.21) 


properties 


(2.22) 


7 


From  Equation  2.20  it  can  be  seen  that  the  <J>  can  thus  be 


derived  if  H(o)  is  known  since 


*(«)  *  JlH{2-Pu>) 

p-i 


Knowledge  of  4>  can  subsequently  be  used  to  find  i|f 


ilro(x)  e  Wq  ,  it  can  be  written  as 


4^0 


Let 


g(n)  =  j‘i|ro(x)<J)o(2x-n) 

-mm 

=  2  ,<t>i(x— |)) 


and  by  a  similar  approach  to  that  used  previously  we 

4ro(x)  = 

13— •  ^ 

m 

=  2  g{n)^o(2x-n) 


(2.23) 


.  Since 


(2.24) 


(2.25) 


can  find 


(2.26) 
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The  Fourier  transform  of  Equation  2.26  is 


T  (2w)  =  GCw)  ®  (w) 


(2.27) 


The  following  properties  result: 

|G(0)|  =  0 

K?(o))|2  +  |G(o+it)|2  =  1  (2.28) 

lf(«)  G(a))  +  H{ta+n)  G(o>+n)  =  0 


Equations  2.22  and  2.28  are  recognized  as  characteristic  of 
"conjugate"  filters  and  are  necessary  to  ensure  orthogonality 

among  the  <j)  and  i|r  families  of  functions.  Equation  2.28 

ensures  that  each  function  of  the  <j>  (commonly  referred  to 

as  the  scaling  function)  is  orthogonal  to  each  function  of  the 

family.  This  last  condition  results  from  the  fact  that 
X  .  [Ref.  4;  pp.  16-23] 

An  example  of  a  function  that  fulfills  the  stated 
conditions  is 

G(«)  =  e'^“if(«+7i)  (2.29) 
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from  which  we  can  find 


gr(n)  =  (-l)^‘”h(l-n)  (2.30) 

Equation  2.30  is  the  defining  relationship  for  describing 
wavelet-scaling  function  dependency  in  the  application 
algorithms.  [Ref.  3:  p.  679] 

2 .  Decomposition  and  Reconstruction  Algorithms 

Recalling  that  c„{n)  =  f)  ,  we  can  derive  the 


decomposition  recursion  algorithm: 


=  fix)  dx 


=  j’[v^J^77(ic)<l>„(x-2"<®'^’rj-2‘®ic)  ]  fix)  dx 

=  v^E  hik)  f<t>„ix-2~”'i2n+k) )  fix)  dx 
* 

=  v^J^i3(Jc)  c„i2n+k) 

* 


(2.31) 


Likewise  we  find 


d^iin)  =  gr(k)  c„i2n+k)  (2.32) 

k 
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To  determine  the  reconstruction  recursion  algorithm  recall 


E  ^in-i<^)‘l>{«-i)Jc+E  («.-!)* 

k-—  k»— 


(2.33) 

(2.34) 


and  substitute  terms  to  get 

~  /  2  35) 

k  k 

Using  a  change  of  variables  the  following  can  be  found 

=  V^-h(^-2ic)  (2.36) 

<<!>««' =  v^flr(n-2A:)  (2.37) 

so  that 

cjn)  =  y/^{"^c^^[k)h{n-2k)+'^d„.^[k)g{n-2k)]  (2.38) 

The  equations  of  this  section  demonstrate  the  pyramidal 
structure  of  multiresolution  analysis  and  readily  lend 
themselves  to  digital  signal  processing  applications.  [Ref  4: 
pp. 25-28] 

The  important  conclusion  is  that,  assuming  knowledge 
of  the  g  and  h  vectors,  the  c  and  d  coefficients  at  any  level 
can  be  completely  determined  from  the  c  coefficients  at  the 
next  higher  level;  also,  the  c  coefficients  at  any  level  can 
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be  determined  from  the  c  and  d  coefficients  at  the  adjacent 
lower  level.  Notice  that  calculations  within  this  pyramidal 

algorithm  do  not  require  explicit  use  of  the  <|>  and  ijf 

functions  since  interlevel  relationships  are  defined  entirely 
by  the  g  and  h  vectors.  Obviously,  the  same  results  would  be 
reached  if  dilations  and  translations  of  the  orthogonal  basis 
functions  were  used  directly  at  each  level  in  the 
decomposition  and  reconstruction,  but  would  be  found  at  great 
computational  cost  because  of  the  requirement  to  calculate 
inner  products. 
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III.  COMPUTER  ZMPLEMEMTATIOM  OF  MULTIRESOLUTIOM  ANALYSIS 


A.  INTRODUCTION 

Our  goal  is  to  investigate  the  use  of  multiresolution 
analyses  to  extract  information  from  an  input  signal.  In  this 
case,  a  one-dimensional  data  sequence  will  represent  the 
input.  To  avoid  having  to  numerically  evaluate  any  inner 
products  we  will  equate  the  Cj)(n)  coefficients  with  the  data 

sequence,  thereby  constructing  an  auxiliary  function  f,  with 
f =ZC(,  (n)<|)o  (n) ,  which  clearly  resides  in  Vq.  Since  CQ(n) 
entirely  describes  the  input,  it  represents  the  highest 

resolution  level  possible  and  the  apex  of  the  pyramidal 
algorithm.  The  multiresolution  framework  previously  described 
can  now  be  used  to  decompose  /  (thus  the  data  sequence  C^in)) 
into  lower  resolution,  i.e.  m  <  0,  approximation  coefficients 
c^(n) ,  and  detail  coefficients  d^(n) .  The  reconstruction 
recursive  equations  can  be  used  to  recover  the  original  f  (and 
therefore  the  data  sequence) .  A  different  multiresolution 
analysis  exists  for  each  scaling  function/wavelet  set.  This 
thesis  is  limited  to  the  Haar  wavelet  and  to  Daubechies*  group 
of  compactly  supported  wavelets.  The  h  vectors  for  these 
wavelets  were  published  in  Reference  1  and  are  listed  in 
Appendix  H  for  convenience.  Describing  the  h  vectors  is  the 
common  method  of  defining  scaling  function/wavelet  sets 
because  all  other  desired  information  can  subsequently  be 
derived  from  them. 
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Several  computer  programs  were  written  to  implement 
various  multiresolution  analyses.  The  programs  were  written 
in  MATLAB  to  take  advantage  of  the  transportability  of  m- 
files  between  various  operating  systems.  The  capabilities  of 
these  programs,  listed  in  Appendices  B-H,  will  be  described  in 
general.  Specific  information  can  be  obtained  by  referring  to 
program  comments  in  the  appropriate  appendix. 

The  normal  entry  point  into  the  collection  of  programs  is 
the  calling  program  wavelet. m  (Appendix  B) .  A  brief 
description  of  the  various  possible  options  is  presented  upon 
entry  and  the  user  can  make  a  decision  based  on  the  choices 
provided.  Specifically,  the  user  may  decide  to  use  the  Haar 
wavelet,  his  own  wavelet,  or  one  from  the  Daubechies  group  of 
nine  compactly  supported  wavelets  and  analyze  with  either  a 
single  phase  or  multiphase  approach. 

The  single  phase  approach  implies  the  decomposition  start 
point  is  the  data  start  point,  i.e.  a  "snapshot”  of  the  entire 
data  sequence  is  processed.  The  term  multiphase  refers  to 
starting  a  decomposition  at  every  possible  data 
sequence/wavelet  phase  relationship,  i.e.  a  "sliding  window" 
approach,  to  maximize  the  signal  analysis  capability.  As 
desired,  wavelets  are  phase  sensitive  but  feature  extraction 
may  be  hindered  if  the  "best"  input  phase  relationship  is  not 
used.  The  multiphase  analysis  therefore  contains  multiple 
sets  of  coefficients,  with  each  set  containing  all  the 
information  necessary  for  reconstruction. 
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B.  SINGLE  PHASE  NULTIRESOLUTION  ANALYSIS 


The  second  program,  haarwave.m  (Appendix  C)  ,  does  the 
"snapshot"  analysis  of  an  input  data  sequence  using  the  Haar 
basis.  Because  of  the  simple  nature  of  the  Haar,  it  is  easy 
to  generate  the  actual  approximation  and  detail  function  at 
each  resolution  level.  The  input  data  is  viewed  as  a 
piecewise  constant  function  (equating  to  the  sample  and  hold 
operation)  to  allow  easy  inner  product  calculation. 
Representations  of  the  properly  weighted  scaling  function  and 
wavelet  at  each  level  are  provided  to  clearly  indicate  the 
dyadic  relationship  between  levels.  Additionally,  the 
reconstruction  of  the  original  sequence  can  be  viewed  as  well 
as  a  plot  of  the  reconstruction  error.  Finally,  distribution 
of  the  energy  of  each  approximation  and  detail  coefficient  and 
the  total  energy  of  each  resolution  level  are  displayed. 
These  energy  distributions  are  the  foundation  for  signal 
analysis. 

The  third  program,  daubwave.m  (Appendix  D) ,  allows  the 
user  to  pick  one  of  the  wavelets  from  Daubechies  group  or  to 
enter  a  valid  user-defined  h  vector  to  define  the  basis  set. 
Basically,  the  same  output  plots  that  were  generated  for 
wvhaar.m  can  be  seen.  Because  of  the  complexity  and 
irregularity  of  the  basis  functions,  however,  only  the  value 
of  the  coefficients  which  are  generated  for  specific  points 
will  be  plotted  instead  of  the  inner  product  representation. 
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C.  MULTIPHASE  WAVELET  ANALYSIS 


The  last  analytical  program,  multiphs.m  (Appendix  E) , 
allows  the  use  of  any  of  the  previously  defined  wavelets  with 
the  multiphase  algorithms.  The  different  coefficient  sets  are 
displayed  simultaneously,  which  is  possible  because  their 
coefficient  indices  interleave  without  interference.  Clearly, 
redundant  information  is  provided  in  the  sense  that  there  are 
more  coefficients  than  would  be  necessary  for  reconstruction 
of  the  data  sequence  over  its  interval  of  support.  However, 
because  of  the  different  zero  padding  necessary  for  each  set 
of  coefficients  a  different  P_f  in  V_  is  defined  for  each  set. 

fn  ffl 

The  user  will  therefore  be  able  to  see  the  P„f  with  the  most 

m 

distinctive  set  of  coefficients  to  ease  analysis. 

D.  SUPPORTING  PROGRAMS 

The  program  basisplt.m  (Appendix  F)  supports  the 
daubwave.m  and  multiphs.m  programs  by  generating  iterative 
plots  of  the  scaling  function  and  wavelet  bases.  It  can  also 
be  called  up  directly  and  will  use  the  vector  "hcoeff”  as  the 

h  coefficients  to  determine  the  4>  and  based  on  a 

graphical  recursion  method. 

The  functions  wvinput.m  (Appendix  G)  and  daubdata.m 
(Appendix  H)  are  called  as  necessary  by  any  of  the 
multiresolution  analysis  programs  to  provide  selected  input 
data  signals  or  the  h  vectors  from  Daubechies*  group, 
respectively.  The  input  signal  options  consist  of  a  sinusoid. 
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a  pseudo-noise  (PN)  sequence,  a  sinusoid  modulated  with  a  PN 
sequence,  and  any  of  the  above  corrupted  by  noise.  Various 
parameters  can  be  modified  to  satisfy  the  user's  needs. 
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IV.  SIGNAL  ANALYSIS 


All  plots  associated  with  this  section  can  be  found  in 
Appendix  A. 

A.  SINGLE  PHASE  ANALYSIS  OF  A  SINUSOIDAL  WAVEFORM 

The  first  signal  to  be  analyzed  will  be  the  sinusoid  shown 
in  Figure  1,  along  with  its  level  0  approximation.  Figures  2- 
5  show  the  Haar  multiresolution  analysis,  stopping  at  level  -4 
since  there  is  no  energy  at  any  lower  level  for  this 
decomposition.  The  reconstruction  plots  (not  shown)  are 
indistinguishable  from  the  decomposition  plots  because  the 
maximum  reconstruction  error  is  15  orders  of  magnitude  below 
signal  levels.  The  energy  contained  at  each  resolution  level 
for  both  approximation  and  detail  coefficients  can  be  seen  in 
Figure  6.  Clearly,  the  maximum  energy  change  occurs  in  the 
detail  coefficients  at  level  -4  and  the  energy  in  the 
approximation  coefficients  at  every  level  is  the  sum  of  the 
energy  of  the  detail  and  approximation  coefficients  from  the 
level  below,  as  expected.  Figures  7  and  8  provide  a  clear 
summary  of  exactly  where  the  signal  and  filter  best  match  with 
respect  to  sample  number  and  resolution  level.  In  this 
example,  the  Haar  decomposition  highlights  the  phase  changes 
associated  with  the  switching  between  positive  and  negative 
half-cycles  because  of  the  signal/wavelet  phase  coherence. 
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The  decomposition  is  repeated  with  a  different 
multiresolution  analysis  based  on  the  6-coefficient  (third 
order)  Daubechies  wavelet.  Only  the  coefficient  energy  plots 
are  shown  in  Figures  9-11  for  comparison  with  the  Haar  plots. 
It  can  be  shown  that  the  time  index  of  coefficients  goes 
beyond  the  original  support  length  of  the  data  (potentially 
out  to  sample  257  in  this  example,  although  the  plots  were 
truncated  at  sample  90  for  display  purposes  and  because  very 
little  energy  is  associated  with  samples  beyond  90)  .  Although 
the  coefficients  outside  the  signal  support  range  are 
generally  very  small,  they  are  necessary  for  completeness. 
Computing  these  coefficients  requires  affixing  zeros  to  the 
sequence,  most  significantly  at  the  trailing  edge.  These 
"edge  effects”  are  minimized  at  the  leading  edge  (no 
coefficients  are  generated  for  negative  time  indices)  by 
shifting  the  h  vector  so  that  nonzero  values  may  initially 
exist  only  over  [- (m-2)  , . . . , 0, 1]  and  then  assigning  the 
coefficient  value  to  the  index  of  the  second  term  from  the 
right  (initially  0) .  This  is  intuitively  satisfying  from  a 
causality  viewpoint  and  also  aids  interpretation  by  only 
having  coefficients  outside  the  data  sequence  on  one  side. 
Leading  edge  effects  consist  of  the  influence  of  leading  zeros 
necessary  for  filter  coefficients  with  negative  time  indices. 
Note  that  edge  effects  increase  as  the  resolution  level 
decreases  because  lower  resolution  coefficients  are  influenced 
by  a  larger  number  of  data  points.  The  extent  of  edge  effects 
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is  also  dictated  by  the  number  of  h  coefficients. 
Coefficients  beyond  the  data  sequence  support  will  not  have 
time  indices  greater  than  [  (2""'-l)  (number  of  h  coefficients  - 
1)  ]  for  any  of  the  decomposition  programs  and  often 
significantly  fewer  (e.g.,  if  the  data  sequence  length  is  a 
power  of  2  the  Haar  decomposition  will  not  generate  any  terms 
beyond  the  data  length,  as  can  be  seen  in  the  previous  plots)  . 

The  next  set  of  plots  (Figs.  12-18)  show  the  phase 
sensitivity  of  the  decomposition  when  the  input  sinusoid  phase 
is  shifted  from  0  to  45  degrees.  The  energy  distribution  is 
no  longer  as  concentrated  as  in  the  previous  set  of  plots  and 
interpretation  is  consequently  more  difficult.  Changing  the 
sampling  frequency  relative  to  the  sinusoid  frequency  will 
have  a  similar  effect. 

B.  SINGLE  PHASE  ANALYSIS  OF  A  BINARY  PHASE  SHIFT  KEYED  SIGNAL 

The  Binary  Phase  Shift  Keyed  signal  shown  in  Figure  19 
will  serve  to  highlight  a  deficiency  in  the  single  phase 
approach.  Notice  with  the  Haar  wavelet  in  Figures  20-22  that 
the  decomposition  is  exactly  the  same  as  the  pure  sinusoid, 
i.e.,  because  of  the  phase  alignment  we  cannot  discriminate 
between  the  two  signals'  coefficient  energy  plots.  The 
Daubechies  6-coefficient  wavelet  decomposition  (Figs.  23-25) 
does  appear  significantly  different  than  the  sinusoid 
decomposition  but  there  is  no  clear  indication  of  every  phase 
reversal.  Thus,  the  results  are  ambiguous. 
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C.  MULTIPHASE  ANALYSIS  OF  A  BINARY  PHASE  SHIFT  KEYED  SIGNAL 

Since  these  responses  would  be  unsuitable  for  signal 
analysis  purposes  the  multiphase  approach  (described 
previously)  was  developed.  It  can  be  viewed  as  starting 
decompositions  at  each  of  the  2"'"  phases  of  each  level  so  that 
the  most  distinctive  representation  can  be  used  for  analysis 
regardless  of  initial  signal  phase.  The  multiphase  plots  are 
shown  in  Figures  26  and  27  for  the  Haar  and  in  Figures  28  and 
29  for  the  Daubechies  6-coefficient  filter.  There  is  clear 
indication  in  both  sets  of  plots  of  the  signal's  phase 
inversions. 

D.  MULTIPHASE  ANALYSIS  OF  A  TRANSIENT  SIGNAL 

Finally,  the  transient  sign?’'.  of  Figure  30  was 
investigated.  Use  of  the  "zoom”  feature  in  the  programs  was 
necessary  to  clearly  se‘>  the  response.  Figures  31-34  show  the 
responses  associated  with  the  Haar  and  Daubechies  6- 
coefficient  wavelets  that  have  been  the  standard  throughout 
this  thesis.  Higher  order  Daubechies  (more  regular)  wavelets 
were  also  experimented  with  in  an  attempt  to  best  fit  the 
signal  and  it  can  be  seen  that  the  18-coefficient  (ninth 
order)  wavelet  used  for  Figures  35  and  36  does  a  good  job  of 
clearly  indicating  the  presence  of  the  signal  and  of 
maximizing  the  concentration  of  energy  into  only  a  few  points 
in  a  single  resolution  level.  This  good  "match"  of  the  signal 
with  the  wavelet  filter  could  be  used  for  signal 
classification  purposes  as  well  as  detection  if  an  a  priori 
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selection  of  the  wavelet  filter  based  on  a  similar  known 
signal  had  occurred. 
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V.  CONCLUSIONS 


Wavelet-based  multiresolution  analysis  is  another  tool  for 
the  signal  analyst  and  great  potential  exists  for  its  use  in 
various  applications.  Signals  that  can  only  be  represented  in 
the  frequency  domain  with  large  numbers  of  significant  terms 
provide  special  motivation  for  wavelet-based  decomposition,  as 
can  be  seen  in  the  transient  signal  example  above.  Signal 
detection  and  pattern  recognition  through  the  use  of  "matched" 
wavelets  may  also  prove  useful,  especially  in  areas  where 
analysis  at  different  resolution  levels  (i.e.,  possible  signal 
dilation/contraction)  is  required.  The  extent  of  the  ultimate 
usefulness  and  popularity  of  wavelets  will  become  known  as 
others  gain  knowledge  of  basic  wavelet  characteristics  and 
incorporate  multiresolution  analyses  into  their  research. 
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APPENDIX  A 


SIGNAL  ANALYSIS  PLOTS 
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MATLAB  Plot  of  Input  Signal  Vector 
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Sample  number  "n"  (Time  =  n  *  0.0625  sec) 


Approximation  function  at  resolution  level  -1 


Figure  2.  Haar  Response  at  Resolution  Level  -1 
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Sample  number  "n"  (Time  =  n  *  0.0625  sec) 


Approximation  function  at  resolution  level  -2 


Sample  number  "n"  (Time  =  n  *  0.0625  sec) 


Approximation  function  at  resolution  level  -3 


Sample  number  "n"  (Time  =  n  *  0.0625  sec) 


)625  sec) 


Normalized  Energy  of  Approximation  Function  vs.  Resolution  Level 


Figure  6.  Haar  Resolution  Level  Energy  Distribution 
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Resolution  Level 


Figure  7.  Haar  Approximation  Coefficient  Distribution 


Figure  8.  Haar  Detail  Coefficient  Distribution 
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Normalized  Energy  of  Approximation  Coefficients 


Figure  9.  Daubechies  Resolution  Level  Energy  Distribution 
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Resolution  Level 


Energy  in  Approximation  Coefficients 


Figure  10.  Daubechies  Approximation  Coefficient  Distribution 
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using  Daubechies  Wavelet  of  order  3 


MATLAB  Plot  of  Input  Signal  Vector 
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Sample  number  "n"  (Time  =  n  *  0.0625  sec) 


Normalized  Energy  of  Approximation  Function  vs.  Resolution  Level 


Resolution  Level 


Energy  in  Approximation  Coefficients  using  Haar 


Figure  15.  Haar  Detail 


Distribution 


Normalized  Energy  of  Approximation  Coefficients 


Figure  16.  Daubechies  Resolution  Level  Energy  Distribution 
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Resolution  Level 


Figure  18.  Daubechies  Detail  Coefficient  Distribution 
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using 


MATLAB  Plot  of  Input  Signal  Vector 


Sample  number  "n"  (Time  =  n  *  0.0625  sec) 


Normalized  Energy  of  Approjdmation  Function  vs.  Resolution  Level 


Resolution  Level 


Figure  21.  Haar  Approximation  Coefficient  Distribution 
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Haar  Detail  Coefficient  Distribution 


Resolution  Level 


Energy  in  ’ 


Energy  in  Approximation  Coefficients 


Energy  in  Detail  Coefficients 


Energy  in  Approximation  Coefficients 


Energy 


using  Da 


xio-3  MATLAB  Plot  of 


Energy  in  E>etail 


Energy  in  Approximation  Coefficients 


using  Daubechies  Wavelet  of  order  3 


Energy  in  Approximation  Coefficients 


Energy  in  Detail  Coefficients 


using  Daubcchies  Wavelet  of  order  9 


APPENDIX  B 


CALLING  PROGRAM  FOR  ACCESSING  MULTIRESOLUTION  ANALYSES 

PROGRAMS  AMD  FUNCTIONS 
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%  Calling  program  for  wavelet  multiresolution  analysis.  The  capabilities 
t  and  options  of  the  entire  set  of  programs  are  presented,  in  general  terms, 
%  in  the  comments  below. 

b=[' 

N1  «  [ '  This  program  will  perform  roultiresolution  analysis  on  input  data' 

'  through  the  use  of  various  orthonormal  bases  of  compact  support.  ' 

'  It  is  interactive  in  nature  and  will  allow  you  the  opportunity  ' 

'  to  select  among  several  different  options  to  best  support  your  ' 

'  specific  needs.  Use  the  Return  key  to  control  movement  through  ' 

'  the  program  (after  prompts  and  questions,  to  view  next  plot,  etc.)' 

»  / 

'  Program  options  include:  ' 

>  / 

'  1.  Selection  of  Basis  Functions  ' 

'  a.  Haar  wavelet  ' 

'  b.  Daubechies  group  of  compactly  supported  wavelets  ' 

'  c.  User-input  wavelet  ("h"  coefficients)  ' 

t  > 

'  2.  Generation  of  Scaling  Function/Wavelet  plots  ' 

»  ' 

'  3.  Selection  of  Input  Data  ' 

'  a.  User-defined  data  vector  (with  optional  noise  background) :  ' 

'  1)  Sine  wave  ' 

'  2)  Pseudo-noise  (PN)  sequence  ' 

'  3)  Sine  wave  modulated  by  PN  sequence  ' 

'  b.  User-input  data  vector  ' 

'  <Return>  for  more  ...'); 

N2  =  ['  4.  Decomposition  algorithm  based  on:  ' 

'  a.  Data  "snapshot"  -  Fixed  start  point;  one  decomposition;  ' 

'  single  phase  approach  ' 

'  b.  "Sliding"  Data  window  -  Floating  start  point;  decomposition  ' 

'  for  each  start  point;  multiphase  ' 

'  approach;  discrete  approximation  to  ' 

'  continuous  wavelet  transform  ' 


'  Do  you  want  to  see;  ' 

/  » 

'  1.  Haar  wavelet  data  "snapshot"  decomposition,  ' 

'  (input  data  treated  as  a  piecewise  constant  function,  ' 

'  allowing  inner  products  to  be  calculated  and  displayed)  ' 

»  t 

'  2.  Daubechies  or  user  wavelet  data  "snapshot"  decomposition,  or  ' 

/  ' 

'  3.  "Sliding"  data  window  wavelet  (any  type)  decomposition?  ') 

clc 

disp(b) 

disp(Nl) 

pause 

clc 

disp(b) 
disp(N2) 
disp (b) 

desire  =  input ( 'Answer  1,  2,  or  3.  '); 
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if  desire  ««  1 
haarwave 
elseif  desire 
daubwave 
else 

nultiphs 

end 


APPENDIX  C 


PROGRAM  FOR  SINGLE  PHASE  ANALYSIS  USING  HAAR  WAVELET 
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%  This  program  performs  single-phase  Haar  wavelet  decomposition  on  a 
%  user-supplied  or  program-generated  input  waveform.  Representation  is 
%  made  with  bar  graphs  to  show  actual  inner  product  calculation  (input 
t  data  is  taken  as  a  piecewise  constant  function) .  The  input  data  sequence 
%  will  be  zero-padded  to  accomodate  the  calculation  of  all  possible  non- 
%  zero  coefficients  and  approximation  (“c")  and  detail  (“d")  coefficients 
t  may  be  computed  for  sample  numbers  beyond  the  data  endpoint.  The  user 
%  should  consider  these  "edge  effects"  in  his  analysis. 

% - 

%  Determine  if  user  wants  to  input  an  external  data  vector  or  desires  to 
%  build  one  through  the  program. 

clc 


b=[' 
Q1  = 


['Do  you  wish  to  use:  ' 

'  1.  your  own  MATIAB  formatted  row  vector,  or  ' 

'  2.  a  program  generated  vector  from  the  following  menu?' 

'  -  Sine  wave  ' 

'  -  Pseudo-noise  (PN)  sequence  ' 

'  -  Sine  wave  modulated  by  PN  sequence 


disp(b) 

disp(Ql) 

Pick  =  input ( 'Answer  1  or  2:  '); 
if  Pick  ==  1, 

%  Read  in  user's  input  vector 

N1  =  ['  Note:  If  the  number  of  data  points  is  a  power  of  2,  the' 

'  results  are  easier  to  interpret  because  there  are  not  ' 

'  any  "edge  effects".  ']t 

disp (b) 
disp(Nl) 

Wavedata  =  input('What  is  the  name  of  your  input  vector?  '); 
sampfreq  =  input ('What  was  the  sampling  frequency  (Hz)?  '); 

Tsample  =  1/sampfreq; 

else 

[Wavedata, Tsample)  =  wvinput (Pick) ; 

end 

I - 

%  Decomposition  algorithm  h(0)=0.5,  h(l)=0.5,  g(0)=»0.5,  g(l)=— 0.5 
datlngth  =  length (Wavedata) ; 

numrows  =  ceil (log(datlngth)/log(2) ) ;  %  find  number  of  resolution  levels 

numpts  =  2'numrows;  Newdata  -  zeros ( 1, numpts) ; 

Newdata(l:datlngth)  -  Wavedata; 

d  =  zeros (numrows, numpts) ;  c  =  zeros(numrows+l, numpts) ; 
detail  =  zeros (numrows, numpts) ;  approx  «  zeros(numrows,numpts); 
c(numrows+l, : )  -  Newdata; 
top  «  numpts; 

const  «=  l/sqrt(2);  %  normalization  constant  *  coefficient  magnitude 

ctr  =  numrows; 
while  ctr  >=  l, 

n  =  numrows-ctr+1 ; 
shift  “  2^n; 
shiftl  =  shift/2; 

top  =  ceil (top/2);  %  number  of  points  at  this  level 
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%  get  ”c"  and  "d"  coefficients 


end 


for  k  =  l:top 
jump  =shift*k; 
indexl  =  jump- (shift-1) ; 
index2  =  jump- (shiftl-1) ; 
firsttnn  »  c(ctr+l, indexl)  ; 
if  index2  <*  niuapts, 

secndtm  =  c(ctr+l,index2) ; 
else 

secndtnn  “0;  %  zero  padding 


end 

d(ctr, indexl)  ■«  firsttnn 
c(ctr, indexl)  =  firsttnn 
for  j  =  0: (shift-1) 
if  j  <  shiftl, 

detail (ctr, indexl+ j ) 
else 

detail (ctr , indexl+j ) 
end 

approx ( ctr , indexl+j )  = 
end 
end 

d(ctr,:)  =  const*d(ctr, ; ) ; 
c(ctr,:)  =  const*c (ctr , : ) ; 
normlize  =  sgrt(2*shift) ; 
detail (ctr ,: )  =  detail (ctr, 
approx ( ctr, : )  =  approx (ctr, 
ctr  =  ctr-1; 


-  secndtnn; 

+  secndtnn; 

%  build  matrices  for  display  purposes 
»  d (ctr, indexl) ; 

*  -d (ctr, indexl) ; 
c(ctr, indexl) ; 


%  normalization  constant  for  this  level 
: ) /normlize; 

:  )/normlizc; 


% - 

%  Plot  "c"  and  "d"  coefficients  by  resolution  level 

indxtoO  =  [0.5:numpts-0.5] ; 

plotmin  =  1.2*min(Mavedata) ;  plotmax  =  l.2*max{Wavedata) ; 

if  ( (plotroin“=0) & (plotmax»“0) ) ,  plotmax  «=  0.5;  end 

if  plotmin  >  0,  plotmin  =  0;  end 

if  plotmax  <  0,  plotmax  =0;  end 

v=[0,numpts, plotmin, plotmax] ; 

axis(v) ; 

bar ( indxtoO, c(numrows+l, ; ) ) 

title ( 'Approximation  function  to  Input  Signal  at  resolution  level  O') 
xlabel ([ 'Sample  number  "n"  (Time  =  n  *  ' ,num2str (Tsample) , '  sec)']) 
pause 
clg 

outptdet  =  zeros(l,l);  outptapp  =  zeros(l,l}; 
for  k  =  numrows:-l:l 

level  =  k-numrows-l; 
step  =  (2^ (-level) )/2 ; 
for  j  =  1 : numpts/step 

outptdet (j)  =  detail (k, (j-1) *step+l) ; 

end 

stepa  “  2* step; 

indxtoOa  =  [ stepa/2: stepa : numpts-stepa/2 ] ; 
for  j  =  l:numpts/stepa 

outptapp(j)  =  approx(k, ( j-1) *stepa+l) ; 

end 

axis(v)  ; 
if  k  >  1, 

subpl ot(211),bar(indxto0a, outptapp) 
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else 

indxtoOa  =  [Oinumpts]; 

approxl  «  approx (1 ,:) ;approxl (numpts+l)  «  approx(l,nu3Bpts) ; 
subpl ot ( 2 1 1 ) , plot ( indxtoOa , approxl ) 

end 

title (( 'Approximation  function  at  resolution  level  ' ,num2str( level) ] ) 
xlabel ([ 'Sample  number  "n"  (Number  of  samples  =  ' ,num2str (numpts) , ' ) ' ] ) 
axis (V) ; 

sxibplot  (212),  bar  ( indxtoO ,  outptdet ) 

title (( 'Detail  function  at  resolution  level  ' ,num2str(level) ] ) 
xlabel ([ 'Sample  number  "n"  (Time  «  n  *  ' ,num2str(Tsample) , '  sec)']) 
pause 
clg 

clear  outptapp  outptdet  indxtoO 
indxtoO  “  indxtoOa; 
clear  indxtoOa 

end 

subplot 

% - 

%  Reconstruction  algorithm 

clc 

recmpout  =  zeros(l,l); 
disp(b) 

N2  =  [ '  Level-by-level  Recomposition  can  be  observed  if  desired.  ' 

'  The  Recomposition  algorithm  starts  with  the  lowest  level  ' 

'  Approximation  Function  and  successively  adds  in  the  Detail' 

'  Function  to  obtain  the  next  higher  level  Approximation.  ']; 
disp(N2) 
disp (b) 

c)calgthm  =  input('Do  you  want  to  see  the  Recomposition  (Y  or  N)?  ','s'); 
if  c)calgthm  ==  'Y', 
level  =  -nurorows; 
axis (v)  ; 

pi ot ( indxtoO , approxl ) 

title ([ 'Level  ' ,num2str( level) , '  Recomposition']) 

xlabel ('Sample  number') 

pause 

clg 

recomp  =  approx ( 1 , : ) ; 
for  K  =  Itnumrows 
level  =  level+1; 
recomp  =  recomp+detail (K, : ) ; 
step  =  2'' (-level)  ; 

indxtoO  »  [ step/2: step: numpts-step/2 ] ; 
for  j  =  l:numpts/step; 

recmpout(j)  =  recomp ( (j-1) *step+l)  ; 

end 

axis (V) ; 

bar ( indxtoO , recmpout) 

title ([ 'Level  ' ,num2str( level) , '  Recomposition']) 

xlabel ( 'Sample  number') 

pause 

clg 

end 

bar ( indxtoO , Newdata-recmpout) 
title ( 'Recomposition  Error') 
xlabel ( 'Sample  number') 
pause 
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clg 

end 

% - 

%  Coefficient  energies  are  used  as  a  measure  of  response 

c  =  c. '2 ; 
d  «  d.^2; 

Enrgytot  -  sum{c(numrows+l, : ) ) ; 

Enrgycro  =  zeros ( 1, numrows+1 ) ;  Enrgydro  »■  zeros ( l, nunrows+1) ; 

for  j  •  l:numrows  %  find  energy  in  each  resolution  level 

Enrgycro(j)  =  (sum{c( j , : ) ) )/Enrgytot; 

Enrgydro(j)  ■=  (sum(d(j , : ) )  )/Enrgytot; 

end 

if  Enrgytot  ==  0, 

Enrgycro  ( numr ows+ 1 )  •»  0 ; 
else 

Enrgycro { numrows+1)  =  1; 
end 

Enrgy dr o ( numr ows+ l )  =  0 ; 

Ivl  =  [-(numrows) : 1: 0] ; 

V  »  [-(numrows+1) ,1,0, 1.2) ; 
axis(v) ; 

subplot(211) , bar (Ivl, Enrgycro) 

title ( 'Normalized  Energy  of  Approximation  Function  vs.  Resolution  Level') 
xlabel ( 'Resolution  Level') 
axis(v) ; 

subplot(212) , bar (Ivl, Enrgydro) 

title ( 'Normalized  Energy  of  Detail  Function  vs.  Resolution  level') 

xlabel ('Resolution  Level') 

pause 

subplot 

clc 

% - 

%  Display  mesh  and  contour  plots  of  coefficient  energy  with  optional  "zoom" 

%  capability  to  aid  analysis 

strtsamp  =  1;  endsamp  «  numpts;  strtrow  -  i;  endrow  «  numrows;  zoom  =  i; 
highres  »  -1;  lowres  -  -numrows; 
while  zoom  ==  1, 

rangea  =  [ -highres-1: -lowres ) ;  ranged  «  (-highres: -lowres] ; 

indxtoO  =  [ strtsamp- 1 : endsarop-1 ] ; 

mesh ( c ( strtrow : endrow+1 , strtsamp : endsamp) ) 

title ('Energy  in  Approximation  Coefficients  using  Haar  basis') 
pause 

contour ( c ( strtrow ; endrow+1 , strtsamp ; endsamp) , 10 , indxtoO , rangea ) 
title ( 'Contour  Map  of  Approximation  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"  (Time  «  n  *  ' ,num2str (Tsample) , '  sec)']) 
ylabel ( 'Decomposition  Number  (=  -Resolution  Level)') 
pause 

mesh ( d ( strtrow : endrow , strtsamp : endsamp) ) 

title ( 'Energy  in  Difference  Coefficients  using  Haar  basis') 
pause 

contour (d (strtrow: endrow, strtsamp: endsamp) , 10, indxtoO, ranged) 
title ( 'Contour  Map  of  Difference  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"  (Time  =  n  *  ' ,num2str (Tsample) , '  sec)']) 
ylabel ( 'Decomposition  Number  (»  -Resolution  Level)') 
pause 
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clc 

showBore  «  input( 'Would  you  like  to  "zoob"  in  on  a  section  (Y  or  N)? 
i f  showBore  ' Y ' , 

disp(b) 

disp('  You  will  set  the  display  saBple  nuBber  and  resolution  level  liaits.') 
disp(b) 

disp(['  SaBple  range  0: ' ,nuB2str (nuapts-l) ] ) 

disp(['  Resolution  range  -1: ' ,nuB2str(-nuBrows) ] ) 

revu  •  input ('Need  to  see  the  original  contour  plots  again  (Y  or  N)?  ','s'); 
if  revu  'Y', 

contour (c, 10, [0:nuBpts-l] , [Oinuarows]) 

title ( 'Contour  Hap  of  ApproxiBation  Coefficient  Energy  Distribution') 
xlabel ([ 'SaBple  nuBber  "n"  (Time  •  n  *  ' ,nuB28tr (Tsample) , '  sec)']) 
ylabel(' Decomposition  Number  (•>  -Resolution  Level)') 
pause 

contour(d, 10, [Oinumpts-l] , [Itnumrows] ) 

title ( 'Contour  Map  of  Difference  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"  (Time  «  n  *  ' ,nuB2str (Tsample) , '  sec)']) 
ylabel ( 'Decomposition  Number  (••  -Resolution  Level)') 
pause 
end 

strtsamp  =  input ( 'Sample  start  point?  ')+l; 
endsamp  =  input ( 'Sample  endpoint?  ')+l; 

highres  «  input ( 'Highest  resolution  level  (least  negative)?  '); 
lowres  »  input (' Lowest  resolution  level  (most  negative)?  '); 
endrow  =  numrows+l+highres; 
strtrow  -  numrows+l+lowres; 
else 

zoom  ■»  0; 
end 
clg 
end 
clc 
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APPENDIX  D 


PROGRAM  FOR  SINGLE  PHASE  ANALYSIS  USING 
DAUBECHIES  OR  USER  WAVELET 
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t  This  program  performs  single  phase  decomposition  of  the  input  waveform 
%  with  either  a  Dauhechies  wavelet  of  user-defined  order  or  with  a  wavelet 
%  provided  by  the  user.  The  input  data  seguence  will  be  zero-padded  to 
%  accomodate  the  calculation  of  all  possible  non-zero  coefficients  and 
%  approximation  (“c")  and  detail  ("d")  coefficients  may  be  computed  for 
%  sample  numbers  beyond  the  data  endpoint.  The  user  should  consider  these 
%  "edge  effects"  in  his  analysis. 


% - 

%  Input  the  "h"  coefficients 

clc 

b  -  [ '  '  J  ; 
disp(b) 

D1  =  [ '  Would  you  like  to  use:  ' 

'  1.  your  own  decomposition  coefficients  (use  wvhaar.m  if  you  desire' 

'  to  use  the  Haar  basis  set) ' 

'  2 .  those  associated  with  Daubechies  compactly  supported  wavelets?  ' ] ; 

disp(Dl) 

hpick  =  input ( 'Answer  1  or  2:  '); 
disp(b) 

if  hpick  ==  1, 

disp('  Input  as  a  row  vector  with  sum  of  coefficients  normalized  to  "1".') 
hcoeff  »  input ('What  is  the  name  of  your  decomposition  coefficient  vector?  '); 
choice  =  1; 
else 

choice  ■=  input('What  is  the  desired  wavelet  order  (2-10  are  available)?  ') ; 
[hcoeff]  =  daubdata (choice) ; 
end 

hlength  •*  length(hcoeff)  ; 

if  hlength  ==  0 
disp(b) 

dispi 'ERROR:  Wavelet  order  selected  is  outside  allc.v'able  range.') 
break 
end 

% - 

%  Determine  if  user  wants  to  see  plots  of  basis  functions 

ckplt  «  input('Do  you  want  to  see  the  Scaling  Function/Wavelet  (Y  or  N)?  ','s'); 
if  ckplt  ««  'Y' 
basisplt 
end 

% - 

%  Get  the  input  data  vector 

clc 

disp(b) 

Q1  ■  ['Do  you  wish  to  use:  ' 

'  1.  your  own  MATIAB  formatted  row  vector,  or  ' 

'  2.  a  program  generated  vector  from  the  following  menu?' 

'  -  Sine  wave  ' 

'  -  Pseudo-noise  (PN)  seguence  ' 

'  -  Sine  wave  modulated  by  PN  sequence  ' ) ; 
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disp(Ql) 

Pick  <•  input ( 'Answer  1  or  2; 
if  Pick  —  1, 

%  Read  in  user's  input  vector 

Wavedata  ■  input ('What  is  the  name  of  your  input  vector?  '); 
sampfreq  «  input('What  was  the  sampling  frequency  (Hz)?  '); 

Tsample  »  1/sampfreq; 
else 

[Wavedata, Tsample]  «  wvinput (Pick) ; 
end 

t - 

%  "g''  coefficients  are  derived  from  the  ”h“  coefficients 

hcoeff  «  sqrt(2)*hcoeff ; 
gcoeff  =  fliplr (hcoeff) ; 
for  j  =  2:2:hlength 

gcoeff (j)  -  -gcoeff (j); 

end 

% - 

%  Decomposition  algorithm 

datlngth  *  length (Wavedata) ; 

numrows  =  ceil (log(datlngth)/log(2) ) ;  %  find  number  of  resolution  levels 

Lastpts  =  zeros (l,numrows+l) ;  Shift  =  ones(l,numrows+l) ; 

Iiastpts(numrows+l)  =  datlngth;  lastpt  =  datlngth; 

for  k  -  numrows: -1:1,  %  find  number  of  points  required  at  each  level 

evnorodd  =  rem (lastpt, 2) ; 
if  evnorodd  ==  0, 

lastpt  *  lastpt/2+ceil( (hlength-2)/2) ; 
else 

lastpt  =  (lastpt-l)/2+fix(hlength/2) ; 
end 

Lastpts (k)  »  lastpt; 

Shift(k)  =  2'' (numrows -k+1)  ; 
end 

hhalf  =  ceil (hlength/2) ; 

dprime  =  zeros(numrows,Lastpts(numrows)+hhalf) ; 
cprime  =  zeros(numrows,Lastpts(numrows)+hhalf) ; 

clastvct  =  zeros(l,datlngth+2*hlength-3) ; 
clastvct(hlength-l:datlngth+hlength-2)  =  Wavedata; 
ctr  “  numrows; 
while  ctr  >  0, 

lastpt  »  liastpts  (ctr)  ; 
nwcvctr  =  zeros(l,lastpt) ; 
nwdvctr  «  zeros(l, lastpt) ; 

for  k  =  1: lastpt  %  coefficients  calculated,  by  resolution  level, 

startpt  «  2*k-l;  %  using  convolution  operation  with  shifts  of  2 

endpt  »  2*k+hlength-2 ;  %  (downsampling) 

nwcvctr(k)  »  hcoef f*clastvct (startpt: endpt) ' ; 
nwdvctr (k)  «  gcoeff*clastvct (startpt; endpt) ' ; 

end 

clastvct  =  [ zeros (l,hlength-2)  nwcvctr  zeros ( 1, hlength-1) ) ; 
cprime (ctr, 1: lastpt)  -  nwcvctr; 
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dpriine(ctr,l;lastpt)  «  nwdvctr; 
ctr  «  ctr-l; 
end 


% - 

t  Coefficient  energies  are  used  as  a  measure  of  response 


cenergy  «  cprime.''2; 
denergy  =  dprime.''2; 


Nl 


The  next  plot  will  show  the  normalized  energy  distributions  of  the' 
approximation  and  detail  coefficients  as  a  function  of  resolution  ' 
level.  Based  on  this  information,  you  will  select  the  minimum  ' 

resolution  level  to  be  displayed  on  all  future  plots.  ' 

/ 

"Edge  effects"  can  cause  coefficients  to  be  generated  for  sample  ' 
numbers  beyond  the  endpoints  of  your  data  sequence  (the  lower  the  ' 
resolution  level,  the  greater  the  sample  number  required;  therefore,  ' 

zero  padding  of  the  original  data  sequence  is  required) .  ' 

/ 

This  program  translates  the  indices  of  the  "h"  and  "g"  vectors  to  ' 
[  -(length  of  vector-2), 1  ].  As  a  result,  coefficients  will  exist  ' 

for  sample  numbers  beyond  the  last  data  point  (n^N) ,  but  will  not  ' 

exist  for  sample  numbers  prior  to  the  first  data  point  (n=0) .  ' 

These  coefficients  are  necessary  for  completeness,  but  may  be  ' 

difficult  to  interpret  because  the  number  of  original  data  points  ' 

used  in  their  calculation  can  be  a  small  percentage  of  the  total  ' 

considered  (the  rest  are  "0").  Similarly,  coefficients  calculated  ' 
for  sample  numbers  near  the  beginning  of  the  sequence  will  also  be  ' 
affected  by  leading  zeros  affixed  to  the  data  sequence.  ' 

9 

Limiting  the  sample  numbers  required  for  display  may  assist  inter-' 
pretation;  select  the  lowest  resolution  level  containing  significant  ' 
energy  for  this  effect.  <Return> 


clc 

disp(b) 

disp(Nl) 

pause 


originlE  =  zeros (1, datlngth) ;  originlE  =  Wavedata.^2; 

Enrgytot  =  sum (originlE) ; 

Enrgycro  =  zeros ( 1, numrows+1 ) ;  Enrgydro  =  zeros ( 1, numrows+1) ; 
c)cerr  =  zeros  ( 1,  numrows)  ; 

for  k  1; numrows  %  find  energy  in  each  resolution  level 

Enrgycro (k)  =  sum (cenergy (k, :)) ; 

Enrgydro (k)  »  sura (denergy (k, :)) ; 

end 

if  Enrgytot  ==  0, 

Enrgycro  (numrows+1)  =»  0; 
else 

Enrgycro (numrows+1)  =  Enrgytot; 

Enrgycro  =  Enrgycro/Enrgytot;  Enrgydro  «»  Enrgydro/Enrgytot ;  %  normalize 

end 

Enrgydro (numrows+1)  -  0; 


ckenergy  »  Enrgycro(l:numrows)+Enrgydro(l:numrows) ; 

t  Energy  balance  check.  The  total  energy  at  any  level  should  equal  the 
t  energy  in  the  "c"  coefficients  in  the  next  higher  level, 
for  k  -  1: numrows 
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ckerr(k)  «  abs(Enrgycro(k+l)-ckenergy(K)) ; 
if  ckenergy(k)  >  10^ (-6) *Enrgytot, 
errenrgy  ••  ckerr(k)/iiiax(Enrgycro(k+l)  .ckenergy  (k) )  ; 
if  errenrgy  >  10^ (-6), 
clc 

disp(b) 

disp('  ERROR:  Energy  check  between  levels  does  not  balance.') 
disp('  If  you  entered  your  own  "h"  vector,  recheck  its  validity.') 

break 
end 
end 
end 

Ivl  =  [-(nunrows) :l:0] ; 
v=  [-(numrows+i) ,i,0,1.2] ; 
axis(v) ; 

subplot(211) ,bar(lvl,Enrgycro) 

title ('Normalized  Energy  of  Approximation  Coefficients') 
xlabel ( 'Resolution  Level') 
axis(v)  ; 

subpl ot(212),bar(lvl, Enrgydr o ) 

title ( 'Normalized  Energy  of  Detail  Coefficients') 

xlabel ( 'Resolution  Level') 

subplot 

pause 

% - 

%  Output  the  number  of  points  required  to  properly  display  coefficients  at 
%  each  resolution  level  (determined  by  "edge  effects") 

clc 

disp(b) 

N2  ■=  ('  Listed  below  are  the  number  of  samples  required  to  properly  display' 

'  each  Resolution  Level .  ' ] ; 

disp(N2) 
disp(b) 

disp('  Resolution  Level  Number  of  samples') 
for  )c  =  isnumrows 

vctrindx  =  numrows-k+1 ; 

nurasamps  =  (Lastpts (vctrindx) -1) ‘Shift (vctrindx) +1; 

disp(['  ' ,num2str{-k) , '  ' ,num2str(numsamps) ] ) 

end 

nmbrlvl  =  -input ('What  is  the  lowest  Resolution  Level  you  desire  to  see?  ') ; 

% - 

%  Plot  "c”  and  "d"  coefficients  by  resolution  level 

indxtoO  “  [0:datlngth-l] ; 

plotmin  =  1.2*min(Wavedata) ;  plotmax  -  1.2‘max(Wavedata) ; 
if  ( (plotmin"0)  i  {plotmax“0) ) ,  plotmax  ■  0.5;  end 
if  plotmin  >  0,  plotmin  -  0;  end 
if  plotmax  <  0,  plotmax  -  0;  end 
V  »  [0, datlngth-1, plotmin, plotmax] ; 
axis(v) ; 

plot ( indxtoO, Wavedata, '*') 

title ( 'Approximation  Coefficients  at  resolution  O') 

xlabel ([ 'Sample  number  "n"  (Time  »  n  ‘  ' ,num2str (Tsample) , '  sec)']) 
pause 

lastrow  «•  numrows-nmbrlvl+1; 
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clmnsize  «  (LAStpts (lastrow) -1) *shift (lastrow)+l; 

c  B  zeros  (runbrlvl-4-1,  clmnslze)  ; 

c(ninbrlvl+l,  Irdatlngth)  •«  originlE; 

d  B  zeros (rmbrlvl, clmnsize) ; 

ctr  «  nunrows;  Ctrl  ■  nmbrlvl; 

while  ctr  >-  lastrow, 

sanpval  •  zeros(l,Lastpts(ctr} } ; 
for  k  =  l:Lastpts(ctr) 

index  *  Shift (ctr) *(k-l)+l; 
c(ctrl, index)  «  cenergy (ctr,k) ; 
d{ctrl, index)  -  denergy(ctr,k) ; 
sanpval (k)  «  index-1; 

end 

plotnin  =  1.2*nin(nin(cprine(ctr, :) ) ,nin(dprine(ctr, : ) ) ) ; 
plotnax  B  1.2*inax(nax(cprine(ctr, : ) ) ,nax(dprime(ctr, : ) ) ) ; 
if  ( (plotnin==0) S (plotnax==0) ) ,  plotnax  =  0.5;  end 
if  plotnin  >  0,  plotnin  =  o;  end 
if  plotnax  <  0,  plotnax  »  0;  end 
V  B  [0,nax(sanpval) .plotnin, plotnax] ; 
n  “  nunrows-ctr+1; 
axis (V) ; 

subplot(211) , plot ( sanpval, cprine (ctr, ItLastpts (ctr) ),'*') 

title ([ 'Approxination  Coefficients  at  resolution  level  ' ,nun2str (-n) ] ) 

xlabel ([ 'Sample  number  "n”  (Last  sample  =  ', num2str (max (sampval) ),')'] ) 

subplot(212) ,plot(sanpval,dprime(ctr,l;I,astpts(ctr) ),'*') 

title ([ 'Detail  Coefficients  at  resolution  level  ' ,nun2str (-n) ) ) 

xlabel ([ 'Sample  number  "n"  (Time  =  n  *  ' ,num2str (Tsample) , '  sec)']) 

pause 

clg 

ctr  =  ctr-1;  Ctrl  =  ctrl-1; 

end 

subplot 

% - 

%  Reconstruction  algorithm 

N3  B  [ '  Recomposition  of  the  original  data  sequence  can  be  checked  if  ' 

'  desired.  The  recomposition  algorithm  starts  with  the  approximation' 
coefficients  at  the  lowest  resolution  level  previously  selected  ' 

'  and  successively  "adds"  in  the  detail  coefficients  of  each  level  ' 

'  until  the  original  resolution  of  the  input  sequence  is  obtained.  ']; 

clc 

disp(b) 

disp(N3) 

ckalgthn  b  input('Do  you  want  to  see  the  Recomposition  (Y  or  N)?  ','s'); 
if  ckalgthn  b*  'y', 

hrecomp  b  fliplr (hcoeff) ;  grecomp  =  fliplr(gcoeff) ;  %  invert  in  time 

cstart  B  cprime (lastrow, :) ; 
ctr  B  lastrow; 
while  ctr  <*  nunrows; 
lastpt  B  Lastpts(ctr) ; 
ckvctr  =  zeros (2, lastpt) ; 

for  k  B  1: fix(hlength/2)  %  compute  higher  level  "c"  coefficients 

ckvctr(l,:)  =■  ckvctr (1, :) +hrecomp(2*k) *cstart (k: lastpt+k-1) +. . 

grecomp(2*k) *dprime(ctr,k; lastpt+k-1) ; 
ckvctr (2 , : )  *  c)cvctr (2, : )+hrecomp(2*k-l) *cstart(k: lastpt+k-1) +. . 

grecomp ( 2 *k-l) *dprime(ctr,k:lastpt+k-l) ; 

end 

if  ren(hlength,2)  o. 
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ckvctr(2 , : ) =ckvctr (2 , ; ) +hrecoinp(hlength) *cstart (hhalf : lastpt+hhalf-1) +. . 
grecoBip(hlength)  *dprinie(hhalf : lastpt+hhalf-1)  ; 

end 

ckvctr  =  reshape(c)cvctr ,  1, 2*Lastpts(ctr) )  ; 
cstart  -  ckvctr; 

lastsamp  =  (Lastpts (ctr+1) -1) ‘Shift (ctr+1) ; 
indxtoO  =  f 0:Shift(ctr+l) ilastsamp] ; 

plotmin  «  1. 2*nin(ckvctr) ;  plotmax  «  1. 2 ‘max (ckvctr) ; 
if  ( (plotniin==0)  S  (plotinax==0) )  ,  plotmax  •=  0.5;  end 
if  plotmin  >  0,  plotmin  «  0;  end 
if  plotmax  <  0,  plotmax  =0;  end 
V  =  [0,  lastsamp, plotmin, plotmax] ; 
axis(v) ; 

plot(indxtoO,ckvctr(l:Lastpts(ctr+l) ) , '*') ; 

title ([ 'Level  ' ,num2str (ctr-numrows) , '  Recomposition']) 

xlabel( [ 'Sample  number  "n"  (Last  sample  =  ' ,num2str (lastsamp) ,')'] ) 

pause 

clg 

ctr  =  ctr+1; 
end 
axis; 

plot ( indxtoO , ckvctr ( 1 : datlngth) -Wavedata , ' * ' ) 
title ( 'Recomposition  Error') 

xlabel ([ 'Sample  number  "n"  (Last  sample  =  ' ,num2str (lastsamp) ,')'] ) 
pause 
clg 
else 
axis; 
end 

% - 

s  Display  mesh  and  contour  plots  of  coefficient  energy  with  optional  "zoom" 

%  capability  to  aid  analysis 

strtsamp  =  1;  endsarop  =  clmnsize;  strtrow  =  1;  endrow  =  nmbrlvl ;  zoom  =  1; 
highres  =  -1;  lowres  =  -nmbrlvl; 
while  zoom  ==  1, 

rangea  =  [ -highres-l: -lowres ] ;  ranged  =  [-highres: -lowres] ; 

indxtoO  =  [strtsamp-l:endsamp-l] ; 

mesh (c ( strtrow : endrow+1 , strtsamp : endsamp) ) 

title( 'Energy  in  Approximation  Coefficients') 

if  choice  -=  1, 

xlabel ([ 'using  Daubechies  Wavelet  of  order  ' ,num2str( choice) ] ) 
end 
pause 

contour ( c ( strtrow : endrow+ 1 , strtsamp : endsamp) , 10 , indxtoO , rangea ) 
title ( 'Contour  Map  of  Approximation  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"  (Time  =  n  *  ' ,nuro2str(Tsample) , '  sec)']) 
ylabel ( 'Decomposition  Number  (=  -Resolution  Level)') 
pause 

mesh ( d ( strtrow : endrow , strtsamp : endsamp ) ) 
title ( 'Energy  in  Detail  Coefficients') 
if  choice  -=  1, 

xlabel ([ 'using  Daubechies  Wavelet  of  order  ' ,num2str (choice) ] ) 
end 
pause 

contour (d (strtrow: endrow, strtsamp: endsamp) , 10, indxtoO, ranged) 
title ('Contour  Map  of  Detail  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"  (Time  =  n  *  ' ,num2str (Tsample) , '  sec)']) 
ylabel ( 'Decomposition  Number  («  -Resolution  Level)') 
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pause 

clc 

showmore  =  input('Would  you  like  to  *'zoob"  in  on  a  section  (Y  or  N)? 
if  showmore  »=  'Y', 
disp(b) 

disp('  You  will  set  the  display  sample  number  and  resolution  level  limits.') 
disp(b) 

disp(['  Sample  range  0: '  ,num2str(cliimsize-l) ) ) 

disp(['  Resolution  range  -1: ' ,num2str(-nmbrlvl) ] ) 

revu  =  input('Need  to  see  the  original  contour  plots  again  (Y  or  N)?  ','s'); 
if  revu  ==  'Y', 

contour(c,l0, [0:clmnsize-l) , [0:nmbrlvl]) 

title ( 'Contour  Map  of  Approximation  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"  (Time  =  n  *  ' ,nuB2str (Tsample) , '  sec)']) 
ylabel ( 'Decomposition  Number  (=  -Resolution  Level)') 
pause 

contour (d, 10, [0:clmnsize-l] , [l:nmbrlvl]) 

title (' Contour  Map  of  Detail  Coefficient  Energy  Distribution') 
xlabel  ([ 'Sample  number  *'n"  (Time  "on*  ',  nuffi2str  (Tsample) , '  sec)']) 
ylabel ( 'Decomposition  Number  (=  -Resolution  Level)') 
pause 
end 

strtsamp  =  input ( 'Sample  start  point?  ')+l; 
endsamp  =  input ( 'Sample  endpoint?  ')+l; 

highres  =  input ( 'Highest  resolution  level  (least  negative)?  '); 
lowres  =  input( 'Lowest  resolution  level  (most  negative)?  '); 
endrow  =  nmbrlvl+l+highres ; 
strtrow  =  nmbrlvl+1+lowres ; 
else 

zoom  =  0 ; 
end 
clg 
end 
clc 
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APPENDIX  E 


PROGRAM  FOR  MULTIPHASE  ANALYSIS  WITH  ANY  WAVELET 
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t  This  program  will  decompose  the  input  waveform  with  a  Haar  wavelet,  a 
t  Dauhechies  compactly  supported  wavelet  of  user-defined  order,  or  with  a 
%  wavelet  provided  by  the  user. 

%  A  "sliding”  data  window  is  used  to  maximize  signal  analysis  features 
%  of  the  decomposition  by  negating  the  effects  of  random  signal  time  of 
%  arrival.  Approximation  and  detail  coefficients  may  therefore  be 
%  computed  at  all  resolution  levels  for  every  "n"  (sample  niunber) . 

%  Because  of  finite  input  data  length,  however,  some  of  the  coefficients 
t  will  suffer  from  "edge  effects",  i.e.,  the  sample  values  beyond  the 
%  end  of  the  data  sequence  are  treated  as  "0"s. 

%  Output  graphs  consist  of  "3-D"  MATLAB  mesh  plots  and  "waterfall" 

%  contour  maps  of  the  approximation  and  detail  coefficients,  in  addition 
%  to  plots  of  the  coefficients  at  each  level  (if  desired) . 

% - 

%  Input  the  "h"  coefficients 

clc 

b  =  [ '  '  ]  ; 
disp(b) 

D1  =  [ '  Would  you  like  to  use:  ' 

'  1.  your  own  set  of  decomposition  coefficients,  or  ' 

'  2.  those  associated  with  the  Haar  wavelet,  or  ' 

'  3.  those  derived  from  Daubechies  group  of  wavelets?']; 

disp(Dl) 

hpick  =  input ( 'Answer  1,  2,  or  3:  '); 

disp(b) 

disp(b) 

if  hpick  ==  1, 

disp('  Kote;  Sum  of  coefficients  must  be  normalized  to  equal  1.') 
disp(b) 

hcoeff=input( 'What  is  the  name  of  your  decomposition  coefficient  vector?  ' ) ; 
choice  =  1; 
elseif  hpick  ==  2, 
hcoeff  =  [0.5  0.5]  ; 
choice  =  1; 
else 

choice  =  input ('What  is  the  desired  wavelet  order  (2-10  are  available)?  '); 
[hcoeff]  =  daubdata (choice) ; 
end 

hlength  =  length (hcoeff ) ; 

if  hlength  *=  0 
disp(b) 

dispi'ERROR:  Wavelet  order  selected  is  outside  allowable  range.') 
break 
end 

% - 

%  Determine  if  user  wants  to  see  plots  of  the  basis  functions 

picture=input ( 'Do  you  want  to  see  the  Scaling  Function/Wavelet  (Y  or  N)?  ','s'); 
if  picture  ==  'Y' 
basisplt 
end 

% - 
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%  Get  the  input  data  vector 


clc 

disp(b} 

Q1  =  ['Do  you  wish  to  use:  ' 

'  1.  your  own  MATIAB  formatted  row  vector,  or  ' 

'  2.  a  program  generated  vector  from  the  following  menu?' 

'  -  Sine  wave  ' 

'  -  Pseudo-noise  (PN)  sequence  ' 

'  -  Sine  wave  modulated  by  PN  sequence 

disp(Ql) 

Pick  •  input ( 'Answer  1  or  2;  '): 
if  Pick  ==  1, 

%  Read  in  user's  input  vector 

Havedata  =  input ('What  is  the  name  of  your  input  vector? 
sampfreq  =  input ('What  was  the  sampling  frequency  (Hz)?  '); 

Tsample  =  1/sampfreq; 
else 

[Wavedata, Tsample)  =  wvinput (Pick) ; 
end 

% - 

%  "g"  coefficients  are  derived  from  the  "h"  coefficients 

hcoeff  =  sqrt(2) *hcoeff ; 
gcoeff  =  fliplr (hcoeff) ; 
for  j  -  2:2:hlength 

gcoeff (j)  -  -gcoeff (j)? 

end 

% - 

%  Decomposition  algorithm 

datlngth  ■  length (Wavedata) ; 

numrows  =  ceil (log(datlngth)/log(2) ) ;  %  determine  number  of  resolution  levels 

Lastpts  =  zeros(l,numrows+l) ; 

for  k  =  1: numrows  %  determine  number  of  points  required  for  each  level 

Lastpts (numrows-k+1)  =  datlngth+ (2^ (k) -1) * (hlength-1) ; 

end 

Lastpts (numrows+1)  «  datlngth;  numpts  =  Lastpts (1) ; 
d  “  zeros (numrows, numpts ) ;  c  =•  zeros(numrows+l, numpts) ; 
c (numrows+1, 1: datlngth)  «•  Wavedata;  cworkvct  »  Wavedata; 
ctr  =  numrows; 
while  ctr  >»  1, 
n  “  numrows-ctr; 

shift  =  2'n;  oldendpt  =  Lastpts (ctr+1) ;  newendpt  -  Lastpts (ctr) ; 

dnewrow  -  zeros (1, newendpt ) ;  cnewrow  -  zeros (1, newendpt) ; 

for  k  «  0:hlength-l  %  coefficient  vectors  are  calculated  for  each 

%  resolution  level  by  adding  shifted,  weighted 
\  versions  of  the  next  higher  level  "c"  vector 
%  (same  effect  as  direct  convolution) 

shiftl  =  k*shift; 

cnewrow (l+shiftl:oldendpt+shiftl)  -  cnewrow (l+8hiftl:oldendpt+shiftl)+. . 

hcoeff  (hlength-k)  *cwor)cvct; 

dnewrow (l+shiftl:oldendpt+shiftl)  =  dnewrow(l+shiftl:oldendpt+shiftl)+. . 

gcoeff (hlength-k) *cworkvct ; 

end 

d (ctr, 1: newendpt)  =  dnewrow; 
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c(ctr, linewendpt)  •  cnewrow; 
cwor)cvct  =  cnewrow; 
ctr  «  ctr-1; 
end 

% - 

%  User  determines  output  format 

clc 

disp(b) 

D3  =  ['  Do  you  want  to  see:  ' 

'  1.  separate  plots  of  the  coefficients  at  each  Resolution  Level' 

'  in  addition  to  the  "3-D"  and  contour  plots,  or  ' 

'  2.  only  the  "3-D"  and  contour  plots?  ']; 

disp(D3) 

pickplot  =  input ( 'Answer  1  or  2:  ') ; 

% - 

%  Plot  "c"  and  "d"  coefficients  by  resolution  level,  if  desired. 

if  pickplot  “=  1, 

indxtoO  =  [0;datlngth-l] ; 

plotmin  =  1.2*min(Wavedata) ;  plotmax  =  l.2*max((favedata) ; 
if  ( (plotmin=-0)  i(plotmax“0) )  ,  plotmax  »  0.5;  end 
if  plotmin  >  0,  plotmin  =  0;  end 
if  plotmax  <  0,  plotmax  °  0;  end 
V  =  [ 0,datlngth-l, plotmin, plotmax ] ; 
axis(v) ; 

plot ( indxtoO , Wavedata , ' * ' ) 

title ('Approximation  Coefficients  at  resolution  O') 

xlabel ([ 'Sample  number  "n"  (Time  =  n  *  ' , num2str (Tsample) , '  sec)']) 
pause 

for  k  -  nuarows:-l;l 
n  =  numrows-k+1 ; 
endpt  “  Lastpts(k); 
indxtoO  =  [0;endpt-l]; 

plotmin  =  1. 2*min(min(c(k, : ) )  min(d (k, ; ) ) ) ; 

plotmax  =  1.2*max(max(c(k, : ) )  ,  IT.  x(d(k, : ) ) ) ; 

if  ( (plotmin==0) i (plotmax==0) )  plotmax  =  0.5;  end 

if  plotmin  >  0,  plotmin  »  0;  enc. 

if  plotmax  <  0,  plotmax  »  0;  end 

V  =  [ 0, endpt- 1, plotmin, pi otroax ] ; 

axis(v) ; 

subplot(211) , plot ( indxtoO, c ( k, 1: endpt) , '*') 

title ([ 'Approximation  Coefficients  at  resolution  level  ' ,num2str(-n) ] ) 
xlabel (( 'Sample  number  "n"  (Number  of  points  -  ' ,num2str (endpt) ,')'] ) 
subplot (212) , plot ( indxtoO, d(k, 1: endpt) , '*') 

title ([ 'Detail  Coefficients  at  resolution  level  ' , num2str (-n) ] ) 
xlabel ([ 'Sample  number  "n"  (Time  -  n  *  ' ,num2str (Tsample) , '  sec)']) 
pause 
clg 
end 

subplot 

axis; 

end 

% - 

\  Coefficient  energies  are  used  as  a  measure  of  response.  Display  mesh  and 
%  contour  plots  of  coefficient  energy  with  optimal  "zoom"  capability  to  aid 
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%  signal  analysis 


c  »  c.*2: 
d  »  d.^2; 

strtsamp  «  1;  endsanp  •«  nunpts:  strtrow  «  1;  endrow  «  niunrows;  zoom  «  1; 
highres  =  -1:  lowres  =  -numrows: 
while  zoom  «=  1, 

rangea  »  [ -highres-1: -lowres ] ;  ranged  -  (-highres: -lowres] ; 

indxtoO  *  [ strtsamp- 1 : endsanp-l ] ; 

mesh  (c (strtrow :endrow-*-l,  strtsamp: endsanp) } 

title( 'Energy  in  Approximation  Coefficients') 

if  choice  -=  1, 

xlabel ([ 'using  Daubechies  Navelet  of  order  ' ,nun2str (choice) ] ) 
end 
pause 

contour ( c ( strtrow : endrow+l , strtsamp : endsamp) , 10 , indxtoO , rangea ) 
title ( 'Contour  Map  of  Approximation  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n”  (Time  »  n  *  ' ,nxm2str (Tsample) , '  sec)')) 
ylabel ( 'Decomposition  Number  (-  -  Resolution  Level)') 
pause 

mesh (d ( strtrow : endrow , strtsamp: endsamp) ) 
title ('Energy  in  Detail  Coefficients') 
if  choice  -=  1, 

xlabel ([ 'using  Daubechies  Wavelet  of  order  ' ,num2str( choice) ] ) 

end 

pause 

contour ( d ( strtrow : endrow , strtsamp : endsamp) ,10, indxtoO , ranged) 

title ( 'Contour  Map  of  Detail  Coefficient  Energy  Distribution') 

xlabel (( 'Sample  number  "n"  (Time  «  n  *  ' ,num2str (Tsample) , '  sec)')) 

ylabel ( 'Decomposition  Number  (=  -  Resolution  Level)') 

pause 

clc 

showmore  =  input( 'Would  you  like  to  “zoom"  in  on  a  section  (Y  or  N)?  ','s'); 
if  showmore  ==  'V', 
disp(b) 

disp('  You  will  set  the  display  sample  number  and  resolution  level  limits.') 
disp(b) 

disp(['  Sample  range  0: ' ,num2str(numpts-l) ) ) 

disp(['  Resolution  range  -l: ' ,num2str( -numrows) ) ) 

revu  =  input('Need  to  see  the  original  contour  plots  again  (Y  or  N)?  ','s'): 
if  revu  «*=  'Y', 

contour (c, 10, (0:numpts-l) , t0:numrows]) 

title (' Contour  Map  of  Approximation  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"  (Time  =  n  *  ' ,num2str (Tsample) , '  sec)')) 
ylabel ( 'Decomposition  Number  -  Resolution  Level)') 
pause 

contour (d, 10, (0;numpts-l) , [l:numrows] ) 

title ('Contour  Map  of  Detail  Coefficient  Energy  Distribution') 
xlabel ([ 'Sample  number  "n"  (Time  -  n  *  ' ,num2str (Tsample) , '  sec)')) 
ylabel ( 'Decomposition  Number  («  -  Resolution  Level)') 
pause 
end 

strtsamp  =  input ( 'Sample  start  point?  ')+l; 
endsamp  «  input ( 'Sample  endpoint?  ')+l; 

highres  =  input ( 'Highest  resolution  level  (least  negative)?  '); 
lowres  »  input( 'Lowest  resolution  level  (most  negative)?  '); 
endrow  »  numrows+l+highres; 
strtrow  •=  numrows+l+ lowres; 
else 
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zoom  0; 
end 
clg 
end 
clc 
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APPENDIX  F 


PROGRAM  FOR  GENERATING  ITERATIVE  PLOTS  OF  WAVELETS 
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\  This  program  creates  and  plots  the  desired  iterative  approximation 
%  to  the  scaling  function  and  the  associated  wavelet  as  determined  by 
%  the  input  "h"  vector.  The  construction  is  based  on  the  "graphical” 
%  recursion  method. 


hcoefnew  »  2*hcoeff; 
m  -  length (hcoefnew) ; 

% - 

%  The  "g"  vector  is  determined  from  the  "h"  vector. 

gcoefnew  =  fliplr (hcoefnew) ; 
for  j  -  2:2:m 

gcoefnew(j)  =  -gcoefnew(j) ; 
end 

% - 

t  Recursively  build  basis  functions 

numbrits  =  input ('How  many  iterations  shall  we  run  for  the  approximation?  '); 
size  =  1; 

hpast  =  ones(l,m); 
newsize  =  m; 

for  i  =  1: numbrits 

%  Build  scaling  function  approximation 

hnew  =  zeros(l, newsize) ; 
for  )c  =  o;size-l 

hnew(2*)c+l : 2*k+m)=hnew (2*k+l: 2*k+m)+hpast (k+1) *hcoefnew; 

end 

%  Get  wavelet  from  scaling  function  approximation  and  also  build  time  vector 

wv  =  zeros (1 , 2*newsize) ;wvlet  =  zeros(l, newsize) ;timevctr  =  zeros (1, newsize) ; 
shift  =  newsize/ (m-1) ; 
for  k  =  0:m-l 

shiftl  =  round (k*shift+l) ; 
shift2  =  round (newsize+k*shift) ; 

wv(shiftl:shift2)  =  wv(shiftl:shift2)+gcoefnew(k+l) *hnew; 

end 

for  k  »  1: newsize 

wvlet(k)  -  wv(2*k-l); 
timevctr(k)  =  k-1; 

end 

interval  »  (1/2) 'i; 

timevctr  ••  interval*timevctr ;  %  scale  time  axis 

s  =  interval*newsize; 

%  Plot  basis  functions  and  check  calculations  with  areas  and  inner  products 
if  i  <»  5, 

plot (timevctr , hnew, '+' , timevctr, wvlet, '*') 

xlabel ([ 'Support  «  ' ,nuro2str(s) , '  (Scaling  Function:  ++++,  Wavelet:  ****)']) 
else 
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plot  (timevctr , hnew, ,tinievctr,wvlet, '  — ') 

xlabel([ 'Support  =  ' ,nujn2str(s) , '  (Scaling  Function:  line.  Wavelet;  dash)'}) 
end 

title ([ 'Approximations  for  Iteration  number  ' ,num2str (i) ] ) 
pause 

hpast  »  hnew; 
size  «  newsize; 
newsize  ■>  2*size+m-2; 

phisuro  sum  (hpast)  ; 
wvsum  «  sum(wlet)  ; 
sp(i)  =  phisum*interval ; 
sw(i)  »  wvsum*interval ; 
ip(i)«  hpast*wvlet' ; 

end 

% - 

t  Output  calculation  chec)cs 

disp('  ') 

disp('  Area  under  Area  under 
disp( 'Scaling  Function  Wavelet 
disp('  ') 

for  i  =  l:numbrits 
disp(['  ' ,num2str(sp(i) ) , ' 

end 

disp('  ') 
disp ( ' <Return> ' ) 
pause 

clear  hcoefnew  gcoefnew  hpast  timevctr  wv 


Inner' ) 
Product ' ) 


' ,num2str (sw(i) ) , '  ' ,num2str (ip(i) ) ] ) 


%  find  area  under  scaling  function 
%  find  area  under  wavelet 

%  find  scaling  function/wavelet  inner  product 
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APPENDIX  6 


FUNCTION  FOR  CONSTRUCTION  OF  INPUT  DATA 
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function  [Data,Tsample]  =  wvinput(x) 


%  This  function  supports  the  various  nultiresolution  programs  by  building 
%  the  input  test  vector  desired  by  the  user  from  the  following  choices: 

%  Sine  wave,  Pseudo-noise  sequence,  and  Sine  wave  modulated  by  Pseudo- 
%  noise  sequence  with  optional  additive  Gaussian  noise. 

% - 

%  Determine  desired  input  signal  type. 

b  =  [ '  ' )  ; 

Q2  =  [ '  Select  the  desired  type  of  waveform:  ‘ 

‘  1.  Sine  wave  ' 

'  2.  Pseudo-noise  (PN)  sequence  ' 

'  3.  Sine  wave  modulated  by  PN  sequence']; 

disp(b) 
disp(Q2) 

Make  =  input ( 'Answer  1,  2  or  3:  '); 

% - 

%  Generate  the  desired  sinusoidal  signal. 

if  Make  -=  2, 
clc 

%  Determine  desired  sine  wave  characteristics. 

N2  “  you  will  determine  the  following  Sine  wave  parameters:' 

'  -  frequency  ' 

'  -  phase  ' 

'  -  amplitude  ' 

'  -  number  of  samples  ' 

'  -  sampling  frequency  or  data  length  ']; 

disp (b) 
disp(N2) 

fc  =  input ( 'Sinewave  frequency  ("fc")  in  Hertz?  '); 
theta  “  input ('Phase  (degrees)?  '); 

A  =  input ( 'Amplitude?  '); 

N  ■=  input ( 'Number  of  Samples  ("N",  with  N  =  a  power  of  2)?  '); 

Q3  «  ['  Do  you  wish  to  set:  ' 

'  1.  the  sampling  frequency,  or' 

'  2 .  the  data  length?  '  ] ; 

disp(b) 
disp(Q3) 

S5  =  input ('Answer  1  or  2.  '); 
if  S5  ==  1, 

N3  =  ('  Note:  Sampling  frequency  (fs)  must  be  selected  so  that  ' 

'  fs/fc  >=  2;  Data  length  will  be  set  to  allow  "N"  samples.']; 
disp(b) 
disp(N3) 

fs2fc  -  input ('Desired  sampling-to-signal  frequency  ratio?  '); 
const  =  2*pi/fs2fc; 
fs  =  fs2fc*fc; 

else 

N4  *  [ '  Note:  Sampling  frequency  will  be  chosen  to  give  "N"  samples' 

'  so  observation  time  must  be  <=  '*N"/(2*"fc") .  ']; 

disp(b) 
disp(N4) 


Tobserv  •  input ( 'Desired  observation  time  (seconds)?  '); 
fs  «  N/Tobserv; 
const  “  2*pi*fc/fs; 

Tsample  •  1/fs; 

end 

Tsample  -  1/fs; 

%  Construct  sine  wave 
theta  “  theta*pi/180; 
for  k  »  1:N 

Data(k)  »«  sin  (const*  (k-1) +theta) ; 

end 

Data  s  A* Data; 

if  Make  ==3,  %  used  if  PN  modulation  is  desired 

Sinedata  °  Data; 
disp(b) 
disp(b) 

N4pt5  “  ['  You  will  new  determine  the  pN  sequence  characteristics.']; 
disp(N4pt5) 

end 

end 


%  This  section  generates  the  desired  pseudo-noise  signal  by  building  the 
%  appropriate  feedback  shift  register  (FSR) . 

if  Make  -=  1, 
clc 

%  Determine  desired  PN  sequence  characteristics. 

N5  =■  ( '  Note:  PN  sequence  =  (+1  or  -l,...,2''m  -  1)  ' 

'  Feedback  shift  register  initial  state  is  {1,1,...,!)']; 
disp (b) 
disp(N5) 

N6  «  ['  You  will  determine  the  following  parameters:' 

'  -  number  of  stages  ' 

'  -  number  and  position  of  taps  ' 

'  -  chip  rate  ' 

'  -  time  delay  ' 

'  -  number  of  samples  ' 

'  -  sampling  frequency  or  data  length  ']; 

disp(b) 
disp(N6) 

mstages  »  input ( 'Desired  number  (m)  of  stages  (2-10)?  '); 
numbrtap  »  input (' Number  of  taps?  '); 
for  k  “  1: numbrtap 

Tap(k) “input ([ 'What  is  the  position  of  Tap  number  ' ,nuo2str(k) , '?  ']); 

end 

chiprate  =  input('Chip  rate  (chips/sec)?  '); 

delay  =  input (' Sequence  delay  in  chips  (positive  real  number)?  '); 
if  Make  «=  2, 

%  If  the  PN  signal  alone  is  desired  the  user  must  input  the  number 
%  of  samples  and  sampling  rate:  otherwise,  they  are  dictated  by 
%  the  choices  made  for  the  sinusoid. 

A  =  sqrt ( 2 ) ; 

N  “  input ( 'Number  of  samples?  '); 

Q4  “  ['  Do  you  wish  to  set:  ' 
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'  1.  the  sampling  frequency,  or' 

'  2.  the  data  length?  ']; 

disp(b) 
disp(Q4) 

P6  =  input ( ' Answer  1  or  2 .  ' ) ; 
if  P6  --  1, 

N7  «  [ '  Note:  Sampling  frequency  ("fs")  must  be  selected  so  that  ' 

'  fs/chiprate  >-  2;  data  length  will  be  set  for  "N"  samples.']; 
disp(b} 
disp(N7) 

fs2cr  “  input ( 'Desired  sampling  frequency  to  chip  rate  ratio?  '); 
Tsample  -  l/(chiprate*fs2cr) ; 
else 

N8  =  [ '  Note:  Sampling  frequency  will  be  chosen  to  give  "N”  samples' 

'  so  observation  time  must  be  <«  "H”/(2*Chip  rate).  ']; 

disp (b) 
disp(N8) 

Tobserv  =  input ( 'Desired  observation  time  (seconds)?  '); 
fs  =  N/Tobserv; 
fs2cr  =  fs/chiprate; 

Tsample  =  1/fs; 
end 
else 

disp(b) 

N9  =  [ '  Same  number  of  samples  ( ' ,num2str(N) , ' )  and  sampling  ']; 

NIO  =  ['  frequency  ( ' ,num2str(fs) , '  Hz)  that  were  selected  for']; 

Nil  =  ['  the  sine  wave  will  be  used.  <Retur.i>']; 

disp(N9) 

disp(NlO) 

disp(Nll) 

pause 

fs2cr  =  fs/chiprate; 
end 

%  Generate  base  PN  sequence 
FSR  =  ones (l.mstages) ; 

L  =  2^nstages  -  1; 
cklength  =  N/(fs2cr*L); 

repeat  =  ceil (cklength) ;  %  Determine  how  many  periods  are  necessary 

V  =  zeros (1 , (repeat+1) *L) ; 

for  k  =  1:L  %  build  one  period  of  the  sequence 

V(k)  =  FSR(mstages) ; 
sigma  =  0; 

for  m  =  l:numbrtap  %  modulo-2  tap  adder 

sigma  *  sigma+FSR(Tap (m) ) ; 

end 

FSR(2 :mstages)  =  FSR(l:mstages-l) ; 

FSR(l)  =  rem(sigma,2) ; 

end 

%  Ensure  enough  periods  are  available  for  the  desired  number  of  samples 
for  n  =  1: repeat 

V(n*L+l: (n+1) *L)  -V(1:L); 

end 

delay  =  rem (delay , L) ; 
for  n  =  1:N 

Data(n)  =  V(fix( (n-l)/fs2cr  +  delay)+l); 
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if  Oata(n)  0,  Data(n)  «=  -1;  end  %  make  sequence  bipolar 

end 

end 


%  If  selected,  modulate  the  sine  wave  with  the  PN  sequence. 

if  Make  *=  3 

Data  Sinedata.  *Oata; 
end 


%  Plot  the  constructed  signal  and  add  noise,  if  desired. 
indxtoO  =  [0:N-1]; 

V  =  [ 0,N-l,-1.2*inax (Data)  ,l.2*inax (Data)  ] ; 
axis (V) ; 

plot (indxtoO, Data) , title ('MATLAB  Plot  of  Input  Signal  Vector') 
xlabel ( 'Sample  number' ) ,ylabel ( 'Volts' ) 
pause 
axis  ; 

clc 

disp(b) 

disp('  White  Gaussian  Noise  is  generated  by  the  MATLAB  "random"  function.') 
noise  =  input ('Do  you  want  noise  added  to  your  signal  (Y  or  N)?  ','s'); 
if  noise  ==  ' Y' , 

%  SNR  based  on  "continous"  signal  energy  and  the  Gaussian  noise  variance 

SNRdB  =  input ('What  is  the  desired  signal-to-Noise  Ratio  (dB)?  '); 

SNR  =  10' (SNRdB/10) ; 

Noisepwr  «  ( (A'2)/2)/SNR; 
rand ( 'seed' , 0) ;  rand ( 'normal ') ; 
noisvctr  =  Noisepwr*rand(l,N) ; 

Data  =  Data+noisvctr ; 

V  =  [0,N-l,-1.2*max(Data) ,1.2*max(Data) ) ; 
axis(v)  ; 

plot (indxtoO, Data) , title ('MATLAB  Plot  of  Input  Data  Vector') 
xlabel ( 'Sample  number') ,ylabel ( 'Volts') 
pause 
axis; 
end 


return 
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APPENDIX  H 


FUNCTION  FOR  DAUBECHIES*  H-COEFFICIENT8 
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function [hcoeff]  =  daubdata(x) 


%  This  function  returns  the  proper  "h"  coefficients  for  the  specified  order 
%  Daubechies  wavelet.  The  input  value  of  "x"  is  the  specified  order. 

if  X  “  2 

hcoeff  »  [0.482962913145;0.836516303738;0.224143868042;-0. 129409522551] ; 
elseif  X  "  3 

hcoeff  -  (0.332670552950;0.806891509311;0.459877502118;-0.135011020010; 

-0. 0854  41273882 ;0. 03  5226291882]  ; 
elseif  X  ==  4 

hcoeff  -  [0.230377813309;0.714846570553;0.630880767930;-0.027983769417; 

-0. 1870348 11719  ;  0. 03 084 1381836 ;0.0328830.i  1667; -0.0105974 01785]  ; 
elseif  X  “=  5 

hcoeff  =  [0.160102397974;0.603829269797;0.724308528438;0.138428145901; 

-0.242294887066; -0.032244869585; 0.077571493840; -0.00624 1490213; 
-0.012580751999; 0.003335725285] ; 
elseif  X  ==  6 

hcoeff  =  [0.111540743350;0.494623890398;0.751133908021;0.315250351709; 

-0.226264693965;-0.129766867567;0.097501605587;0.027522865530; 
-0. 03 15820393 18;  0. 000553842201; 0.004777257511  ,--0.0010773 01085]  ; 
elseif  X  ==  7 

hcoeff  =  [0.077852054085;0.396539319482;0.729132090846;0.469782287405; 

-0.143906003929;-0.224036184994;0.071309219267;0.080612609151; 
-0. 03  8029936935  ;-0. 016574  541631  ,-0.012550998556  ,-0.0004  29577973  ; 
-0.001801640704,-0.000353713800]  ; 
elseif  X  ==  8 

hcoeff  =  [0.054415842243,-0.312871590914,-0.675630736297,-0.585354  683654; 

-0.015829105256:-0.284015542962;0.000472484574,-0.128747426620,- 
-0.017369301002;-0.044088253931,-0.013981027917,-0.00874  6094047,- 
-0.004870352993; -0.000391740373 ,-0.00067544  94  06 ,--0.000117476784]  ; 
elseif  X  ==  9 

hcoeff  =  [0.038077947364,-0.243834  674613;0.604823123690,-0.657288078051,- 

0. 13  3 19738  5825  ,--0.293273783279; -0.09684  07832  23  ,-0.148540749338; 
0.030725681479,--0.067632829061;0.000250947115,-0.022361662124; 
-0.004  723204758,--0.004  281503682;0.00184764  6883,-0.000230385764,- 
-0.000251963189,-0.00003  9347320]  ; 
elseif  X  ==  10 

hcoeff  =  [0.026670057901,-0.188176800078,-0.527201188932,-0.6884590394  54  ; 

0.281172  3  4  3661,--0.24984  64  24  327,--0.195946274  377,-0.12736934  03  3  6,- 
0.093057364  604,--0.071394147166;-0.029457536822,-0.033212674  059,- 
0.003606553567,--0.010733175483,-0.001395351747,-0.001992405295,- 
-0.000685856695; -0.0001164  66855 ,-0.000093588670 ,--0.000013264203]  ; 

else 

hcoeff  =  [  ]  ,- 
break 
end 

hcoeff  =  hcoeff '/sqrt (2) ;  %  normalize  the  "h"  vector 

return 
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